Combined ReaxFF and Ab Initio MD Simulations of Brown Coal Oxidation and Coal–Water Interactions

In this manuscript, we use a combination of Car–Parrinello molecular dynamics (CPMD) and ReaxFF reactive molecular dynamics (ReaxFF-MD) simulations to study the brown coal–water interactions and coal oxidation. Our Car–Parrinello molecular dynamics simulation results reveal that hydrogen bonds dominate the water adsorption process, and oxygen-containing functional groups such as carboxyl play an important role in the interaction between brown coal and water. The discrepancy in hydrogen bonds formation between our simulation results by ab initio molecular dynamics (CPMD) and that by ReaxFF-MD indicates that the ReaxFF force field is not capable of accurately describing the diffusive behaviors of water on lignite at low temperatures. The oxidations of brown coal for both fuel rich and fuel lean conditions at various temperatures were investigated using ReaxFF-MD simulations through which the generation rates of major products were obtained. In addition, it was observed that the density decrease significantly enhances the generation of gaseous products due to the entropy gain by reducing system density. Although the ReaxFF-MD simulation of complete coal combustion process is limited to high temperatures, the combined CPMD and ReaxFF-MD simulations allow us to examine the correlation between water adsorption on brown coal and the initial stage of coal oxidation.


Introduction
Coal continues to be the second largest fuel, which supplies 27% of global primary energy needs [1]. The major consumer of coal resources is coal-fired power plants. In addition, brown coal is one of the most widely used combustion fuel sources in electric power generation, which can be attributed to its abundant reserves and low mining cost [2]. However, the high inherent moisture of low-rank coals (e.g., lignite) obstructs the largescale utilization of low-rank coals [3], as the high water content degrades the energy quality of the coal and results in extra costs for transportation. Hence, upgrading brown coal by reducing its water content is necessary [4]. Moreover, brown coal with complex carbonaceous structures is at a high risk of spontaneous combustion, since brown coal has a large amount of oxygen-containing polar functional groups (including carboxyl, hydroxyl, and carbonyl) [2,5]. Therefore, revealing the relationship between brown coal chemical structures and its oxidation at various temperatures can throw some light on brown coal spontaneous combustion and gasification [6]. In addition, those polar functional groups also play an important role in water adsorption [2,3,7]. Thus, it is of great interest to investigate how the water adsorption is affected by the chemical structure of brown coal.
Extensive ab initio calculations (e.g., DFT) have been carried out to investigate the interactions between coal and water molecules [3,8,9]. Gao et al. [8] obtained several typical conformations of water-lignite complexes and corresponding electrostatic potential distributions using DFT calculations. In addition, they have observed that hydrogen bonding dominates the water adsorption on lignite surface, in which water molecules tend to aggregate around oxygen-containing polar functional groups. Xiao et al. [9] examined the interactions between water molecules and oxygen/nitrogen atoms in a lignite model molecule with a DFT and AIM study, and they demonstrated that water clusters play an important role in water-lignite interaction. Although such ab initio studies of brown coal-water interaction provide detailed structural information about water adsorption on lignite, it is still challenging to access the dynamical information about water transportation within brown coal under different conditions. Additionally, the information about entropy is lacking in those calculations. On the other hand, some coarser models based on percolation theory have been developed to study the coal combustion process [10] as well as coal pyrolysis [11], using the coal's industrial analysis value. To bridge the gap between such crude models and the computationally expensive ab initio studies, some molecular dynamics simulation techniques which can be used to explore the reactive systems might be required. More importantly, since the water adsorption on lignite, water/oxygen transportation within lignite, and lignite oxidation are all correlated, reactive molecular dynamics simulations might allow us to investigate the coal desiccation, spontaneous combustion, and gasification without arbitrary assumptions. Furthermore, reactive molecular dynamics simulations have an advantage over the equilibrium modeling of coal oxidation (e.g., gasification [12]), which employs the free energy minimization method in investigating the systems far from equilibrium.
In the last two decades, the ReaxFF reactive force field developed by van Duin et al. [13,14] has wide applications in simulating high energy materials [15][16][17], hydrocarbon oxidation [18], fracture of graphene [19], and thermal fragmentation of fullerenes [20], etc. ReaxFF relates the bond order/bond energy to separations between atoms so that bonds' formation/dissociation can take place at an appropriate distance. In addition, the parameterization of ReaxFF reactive force fields was based on the DFT calculations to determine bond energy [13]. Since the ReaxFF reactive molecular dynamics is about two orders of magnitude faster compared to quantum mechanics calculations [13,21,22], ReaxFF in combination with proper coal models [23,24] has a great potential in simulation of coal oxidation. During the last decade, ReaxFF was used to examine coal pyrolysis [25][26][27][28][29][30][31], lignin pyrolysis [32], petroleum coke pyrolysis [33], coal combustion [22,34], coal gasification [6,35], and coal tar pitch oxidation [36]. The decomposition of coal molecules to produce small fragments [34] as well as pathways of oxidation [35] have been carefully studied by ReaxFF-MD. However, most of those studies focus on high-temperature oxidation/reactions of coal molecules. Hence, the effects of interaction between water/oxygen and coal on low-temperature coal oxidation have not been fully explored by simulations at a molecular level.
Many studies on diffusion of different gases on coal using molecular dynamics simulations with conventional nonreactive force fields (e.g., COMPASS [37,38], Dreiding [39]) have been reported recently. Although the adsorption and diffusive behaviors of water, oxygen, and other gases on coal can be accurately captured through such conventional molecular dynamics simulations, the coupling between adsorption and oxidation cannot be investigated in such MD simulations, since the covalent bonds are unbreakable in those conventional force fields. A recent research on ab initio molecular dynamics simulation of lignite-water interaction has been carried out [40], which might overcome the limitations of nonreactive MD. Using the Car-Parrinello molecular dynamics simulations, Gao et al. [40] explored the adsorptions of water molecules on lignite surface and observed the replacement phenomenon of the water during a very short period of time. Moreover, their research can provide information about statistical properties of the systems which is lacking in previous DFT studies [3,8,9]. In addition, the Car-Parrinello molecular dynamics simulations results can be useful in verifying ReaxFF-MD simulations, since the training dataset for ReaxFF parametrization does not contain all possible reactions.
In this manuscript, a combination of ab initio molecular dynamics (Car-Parrinello molecular dynamics) simulations and ReaxFF reactive molecular dynamics simulations was adopted to investigate brown coal-water interaction as well as coal oxidation at various temperatures. As both DFT studies [8,9] and ab initio molecular dynamics simulations [40] demonstrated that the hydrogen bonds play a very important role in brown coal-water interactions, the hydrogen bond formation between the molecular model of brown coal and water molecule was explored by both Car-Parrinello MD and ReaxFF-MD in this work. In order to examine the brown coal oxidation under different conditions, only ReaxFF-MD simulations were performed, for the Car-Parrinello MD simulations are too computationally expensive for large systems. Although it is of great interest to study how the accumulation of heat during low-temperature coal oxidation facilitates spontaneous combustion of coal, it is still extremely difficult to simulate such process with relatively fast ReaxFF-MD. Because it is necessary to remove the thermostat for ReaxFF-MD simulation so as to model the heat accumulation when low-temperature coal oxidation occurs. As a result, the temperature increase of the system only relies on the accumulation of heat caused by rare events such as low-temperature coal oxidation. Obviously, enough heat accumulation (e.g., 100 degrees temperature increase) is only possible for large systems which can hardly be studied by molecular dynamics simulations. Our strategy in this manuscript is to simulate the brown coal-oxygen/water systems under different conditions. Hopefully, a series of ReaxFF-MD simulations results can provide us with clues to elucidate spontaneous combustion of brown coal. A more accurate description of spontaneous combustion of coal by ReaxFF-MD might require the application of extreme value theory along with ReaxFF-MD simulations.

Simulation Models and Methods
In this work, the representative structure of lignite (as shown in Figure 1) developed by Meng et al. [41] was employed to simulate lignite combustion and gasification at various temperatures using ReaxFF molecular dynamics simulations [15][16][17]. Both the chemical structure and the relaxed configuration of the lignite molecular model are given in Figure 1. Ten brown coal molecules were placed in periodic boxes of various sizes (e.g., the simulation system as depicted in Figure 1c). Different numbers of water molecules and/or oxygen molecules were then added into the system to investigate their effects on brown coal oxidation. Detailed information about our simulation systems is presented in Table 1. In addition, the initial configurations of all our ReaxFF-MD simulations were prepared by randomly positioning different numbers of those molecules given in Table 1 in the simulation box, according to their concentrations. Simulation systems were first equilibrated by performing a 10 ps constant-temperature, constant-volume (NVT) simulation at 300 K. In addition, then, the simulation systems were heated to the target temperatures rapidly, followed by another 250 ps NVT simulation except where otherwise stated. Such simulations were performed to investigate the chemical reactions in our simulation systems with a 0.1 femtosecond ( f s) timestep, the same as in the equilibration stage. All ReaxFF-MD simulations of brown coal oxidation in this work were performed using the REAX package [13,15,16] of LAMMPS [42] as well as ADF software [43]. Visualization of molecules ( Figure 2c) as well as calculation of various coordinate root-mean-squared deviation (RMSD) were carried out by VMD [44].
To compare the simulation results of ReaxFF-MD to that of Car-Parrinello molecular dynamics directly, a small fragment of the original brown coal model molecule (as shown in Figure 2a) was adopted for low-temperature simulations using both ReaxFF-MD and CPMD, as ab initio molecular dynamics (CPMD) is much more computationally expensive than ReaxFF-MD. For all low-temperature simulations, one model molecule (C 17 H 25 O 4 N) with eight water molecules were randomly positioned in a 14.4 × 12.0 × 12.0 Å 3 periodic box, followed by 10 ps NVT ReaxFF-MD equilibration at 300 K. Simulations with this equilibrated initial configuration were performed using both ReaxFF-MD as well as CPMD, as the density of the system was fixed to 0.3565 g/cm 3 . Car-Parrinello molecular dynamics simulations were carried out using DFT calculation with Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [45] and Vanderbilt ultrasoft pseudopotentials [46][47][48][49]. As suggested by Vanderbilt et al. [49], a plane-wave cutoff of 40 Ry and a charge-density cutoff of 200 Ry were used for all ab initio molecular dynamics calculations in this manuscript.
In addition, the NVT CPMD simulations equipped with a Nosé-Hoover thermostat were run with a time-step of 0.12 f s and a fictitious electronic mass of 400.0 a.u. All the CPMD simulations were performed using the Quantum Espresso package [50][51][52].   [41]. The carbon, oxygen, hydrogen, and nitrogen atoms are represented by black, red, white, and blue beads, respectively; (b) equilibrated chemical structure of the model lignite molecule; (c) an example simulation box of brown coal oxidation (S12 in Table 1).

Lignite-Water Interaction at Low Temperatures, CPMD vs. ReaxFF-MD
As discussed above, low-rank coal, such as brown coal, contains high moisture content. To examine the adsorptions of water molecules on brown coal (especially on the oxygencontaining polar functional groups), the interactions between brown coal model molecule and water molecules were simulated using both Car-Parrinello molecular dynamics and ReaxFF-MD. For CPMD, both relaxation of electronic degrees of freedom and structural optimization (ionic and electronic relaxation) were carried out before MD runs. As shown in Figure 3, the DFT total energy reaches a plateau after structural optimization. Then, a 3-ps CPMD simulation was performed to equilibrate the system for different temperatures, which was followed by a 30-ps CPMD simulation at the same temperature. For example, the system temperature is plotted against simulation time for NVT ensemble at 400 K, as demonstrated in Figure 2b. For comparison, a 30-ps ReaxFF-MD with the same initial configuration was run at 300 K. As demonstrated in Figure 4a, the root-mean-squared deviations (RMSD) of all atoms coordinates in simulation box are plotted against a time lag for both ReaxFF-MD as well as CPMD at 300 K. Using the TRAVIS package [53,54], the average diffusion coefficients (D avg ) of all atoms at 300 K for both ReaxFF-MD and CPMD are obtained which are 1.89 × 10 −9 m 2 /s and 1.53 × 10 −9 m 2 /s, respectively. The good agreement between RMSD/D avg obtained by CPMD and that obtained by ReaxFF-MD simulation shows that the average diffusive behavior of all atoms at low temperature is not affected by dynamical simulation methods. Note that the density of our systems for low temperature calculations was fixed to 0.3565 g/cm 3 , which is much smaller compared to the density of lignite (∼1.2 g/cm 3 ). Recent ReaxFF molecular dynamics simulations of coal pyrolysis and combustion [6,22,28,35] and CPMD simulations of coal-water interactions [40] mainly focused on coal models with low densities ranging from 0.1∼0.7 g/cm 3 [6,22,35]. Therefore, our simulation results can be directly compared to those previous results.
By tracking the trajectories of all atoms during CPMD simulations, a hydrogen bond between the carboxyl group and a water molecule was detected. In addition, the distance between the oxygen of the water molecule and the hydrogen of the carboxyl group is plotted against simulation time as shown in Figure 4b. It can be clearly seen from Figure 4b that, once the hydrogen bond formed, the adsorption of water on carboxyl group remained unchanged for the rest of the simulation. Note that the van der Waals radius of hydro-gen atom is 1.06 Å and the van der Waals radius of oxygen atom is 1.42 Å [40,55]. This observation does not agree with the recent simulation results of CPMD with Vanderbilt pseudopotentials obtained by Gao et al. [40], since Gao et al. observed replacement phenomenon of the water during their 30-ps Car-Parrinello molecular dynamics simulations. This disagreement between simulations results might be partly attributed to the fact that Gao et al. [40] chose 25 Ry as the kinetic energy cutoff for wavefunctions instead of 40 Ry, which was suggested by Vanderbilt [48,49], the developer of Vanderbilt pseudopotentials. In addition, Gao et al. [40] have not specified the kinetic energy cutoff for charge density, and the distances between multiple pairs of oxygen and hydrogen exceeded the half-length of the periodic box in their simulations [40], which might be caused by an incorrect analysis.
In order to investigate how the simulation methods and temperature affect hydrogen bond formation in gas phase for our lignite-water system, the average number of hydrogen bonds for both CPMD and ReaxFF-MD are plotted in Figure 4c. In addition, 3.2 Å and 3.0 Å were selected as the statistical distance of hydrogen bonds (H-O· · ·H) to determine whether a hydrogen bond formation took place. As demonstrated in Figure 4c, the average number of hydrogen bonds obtained by ReaxFF-MD is smaller than that obtained by CPMD, which implies that the hydrogen bonds are not accurately described by ReaxFF. In addition, ReaxFF cannot describe the London dispersion (van der Waals (vdW) attraction) accurately [56], which suggests that ReaxFF might not be sufficient to accurately simulate the coal-water interactions at low-temperature. One possible solution to this problem of ReaxFF was proposed by Liu et al. [57], who introduced the low-gradient model to account for the van der Waals attraction. However, their proposed ReaxFF-lg model was based on calculations at room temperature, which might fail at high temperature. Hence, it remains challenging to simulate spontaneous combustion of coal, for an improved unified reactive force field to describe coal oxidation for a wide temperature range is necessary.

Potential Energy vs. Time for ReaxFF-MD Simulations
The influence of time step on ReaxFF-MD simulations of hydrocarbon reactive gases has been systematically studied by Jensen et al. [58]. Their results demonstrated that system potential energy converged when the time step is shorter than or equal to 0.1 f s. Therefore, for this manuscript, 0.1 f s was adopted as the time step for ReaxFF-MD simulations. In addition, the potential energy of simulation system (S2 in Table 1) was plotted against simulation time at various temperatures in Figure 5. As shown in Figure 5, the potential energy decreases as time increases, for the combustion/oxidation of brown coal with oxygen is exothermic [22]. In addition, potential energy decreases more rapidly as temperature increases in the early stage, due to the increase in reaction rate. Because the combustion process takes much time [34], a 250-ps ReaxFF-MD simulation is not sufficient to equilibrate the system. For low temperature (500 K), the simulated system reaches equilibrium state rapidly because coal combustion does not occur at this temperature, which is shown by Figure 6. As the rate of coal combustion is enhanced by temperature increase, the potential energy curves keeps going down after 200 ps for 1000∼2000 K systems, while the potential curves above 2000 K reaches a plateau within 200-ps. In other words, systems at higher temperatures achieve the equilibrium state faster.  Table 1) as a function of simulation time at different temperatures.
Reaction mechanisms of brown coal oxidation at different temperatures were analyzed by ADF software [43] with Döntgen's method [59]. Thousands of chemical species and elementary reactions were observed during our ReaxFF-MD simulations. Several examples of elementary chemical reactions of S2 system in Table 1 at 3000 K are depicted in Figure 7. Similar to the findings of previous studies [6,22,35], our simulation results demonstrated that the brown coal oxidation was initialized by thermal decomposition to produce smaller fragments. Nevertheless, the charge-equilibration (QEq) calculation adopted by ReaxFF-MD has problems that cause unreasonable charge distribution predictions for chemical structures far from equilibrium [56]. As a result, some of our predicted elementary chemical reactions at high temperature might be unreal. Thus, in this manuscript, we focus on CO/CO 2 /H 2 /H 2 O formation which has been extensively investigated using ReaxFF-MD simulations [6,15,16,22,35].  (Table 1); (f-j): S2 system (Table 1).

Product Analysis of ReaxFF-MD Simulation Results-Fuel Rich
As described above, a series of ReaxFF-MD simulations of brown coal oxidation were carried out to investigate how the products of brown coal oxidation were affected by temperature as well as coal density. In addition, such simulations were performed for both fuel rich and fuel lean conditions. As demonstrated in Figure 6, the concentrations of major products of S1 and S2 (as shown in Table 1) such as H 2 O, CO, CO 2 , and H 2 are plotted against simulation time. For S1 & S2 systems, 10 model molecules of brown coal react with 100 oxygen molecules. As shown in Figure 6a,f, both S1 and S2 exhibit sharp decreases in oxygen numbers above 2500 K. On the other hand, the oxygen molecule numbers only decrease slightly when system temperature is less than or equal to 1000 K within 250-ps simulations. Note that, for relatively low temperature (≤1000 K), oxygen molecules are absorbed to brown coal molecules (or fragments) to form some intermediate products, where the final products like CO/CO 2 /H 2 are negligible (as shown in Figure 6c-e,h-j).
As system temperature increases from 500 K to 3000 K, the number of water molecules increases more and more rapidly over time. For both S1 & S2, the number of water molecules reach a plateau within 50 ps, which seems to suggest that chemical equilibrium is achieved. Nevertheless, as discussed above, the systems cannot be equilibrated within 250-ps ReaxFF-MD simulations even for 3000 K simulations, which can be seen from Figure 6c,e,h,j that CO & H 2 concentrations keep on increasing during the simulations. As shown in Figure 6g, the curve that shows the relationship between water molecule number and simulation time bends downward after 200 ps simulation time, which indicates that water molecules and carbonaceous structures react to form hydrogen molecules in this stage. For high-temperatures (≥2500 K), the number of carbon dioxide as products increases rapidly in the initial stage and then decreases gradually after 50 ps for fuel rich condition, since CO 2 is not the preferable product in this condition. The only difference between S1 and S2 is that the density of S2 is only one half of S1 system density. Due to the entropy gain, the number of small inorganic molecules such as CO and H 2 increases significantly as density decreases, as shown in Figure 6. This fact raises the question whether the reaction pathways observed in previous ReaxFF-MD simulations are applicable to other conditions, especially considering that most previous ReaxFF-MD studies focused on low densities ranging from 0.1-0.7 g/cm 3 [6,22,35].
The reaction between brown coal model molecules and water at various temperatures have also been examined by ReaxFF-MD simulation in this work. The numbers of major products of S3 & S4 systems (in Table 1) as functions of simulation time are presented in Figure 8. The number of water molecules almost remains a constant at low temperatures. As temperature increases, the number of water molecules increases over time in the early stage and then enters a plateau region. For the S4 system at 3000 K (in Figure 8e), the number of H 2 O first increases and then decreases over time, which implies the reaction between water molecules and brown coal to produce hydrogen molecules dominates after 50 ps for the fuel rich condition. Due to the lacking of oxygen molecules, CO production is favored compared to CO 2 production, especially at high temperatures. Hence, although the numbers of water molecules and CO 2 molecules vary insignificantly, the numbers of CO and H 2 increase rapidly above 2500 K. Again, due to the entropy gain as a result of increase of simulation box size, the numbers of CO and H 2 molecules increase more quickly as system density is reduced by half. As demonstrated in Figure 9, the concentrations of major products of S5 & S6 systems (in Table 1) at different temperatures obtained by ReaxFF-MD simulations are plotted against simulation time. As oxygen and water molecules coexist in the S5 & S6 systems for the initial configuration, the number of oxygen molecules decrease faster than that of S1 & S2 systems, which suggests that the interactions among brown coal, water, and oxygen might be cooperative. Figure 9b,g demonstrate that the number of water molecules increases rapidly during the first 50 ps simulation and then decreases gradually for high temperatures, which indicates that water molecules react with lignite to produce hydrogen molecules in this condition. Due to insufficient oxygen present in this combustion reaction, incomplete conversion takes place for this fuel rich condition where a large amount of CO molecules have formed as shown in Figure 9c,h. Note that our simulation results (Figure 9d,i) indicate that CO 2 production is favored at lower temperature (2000 K), whereas CO production is favored at higher temperature (3000 K), which agrees with previous research [22]. In addition, hydrogen molecules only form at high temperatures (2500∼3000 K). When the system density is reduced by half, the amount of gaseous products like CO, H 2 increase significantly, which again confirms the important effects of density on brown coal oxidation.  (Table 1); (e-h): S4 system (Table 1).  (Table 1); (f-j): S6 system (Table 1).
For comparison, the reactions between 10 brown coal model molecules and 1000 water molecules (S9 & S10) were also investigated using ReaxFF-MD simulations. In addition, the amounts of gaseous products as functions of simulation time are plotted in Figure 10. Since there are 1000 water molecules in S9 & S10 systems, the number densities of brown coal molecules in S9 & S10 systems are less than that in S3 & S4 systems. Changes in numbers of water molecules in the S9 & S10 systems are relatively small except S10 at 3000 K. For S10 at 3000 K, the number of water molecules decreases significantly over time, which indicates that the reaction between water and brown coal to produce hydrogen molecules is sped up appreciably at 3000 K. By comparing Figure 8b to Figure 10b (or, Figure 8d vs. Figure 10d, Figure 8f vs. Figure 10f, Figure 8f vs. Figure 10f), we can conclude that the production rates of CO/H 2 molecules in S9 & S10 systems are similar to that in S3 & S4 systems. For this fuel rich condition, the CO 2 production can be neglected compared to CO production. Moreover, by reducing the system density by half, the number of hydrogen molecules increase from ∼140 to ∼200 within the 250 ps ReaxFF-MD simulations.  (Table 1); (e-h): S8 system (Table 1).

Product Analysis of ReaxFF-MD Simulation Results-Fuel Lean
The amounts of gaseous products as functions of simulation time are plotted in Figures 11 and 12 for fuel lean conditions. Due to the larger number of O 2 molecules present in systems, the rate of product generation is greater in comparison to fuel rich condition, which is evident by Figure 11d,h. Such observation in our simulations also agrees with previous experimental results [60,61] and simulation results [22]. The excess of oxygen molecules in S7 & S8 allows the complete combustion of brown coal to take place. As shown in Figure 11d,h, the number of CO 2 increases rapidly at the beginning and then enters a plateau region within 250 ps for high temperatures (2000∼3000 K), since most carbon atoms are changed into carbon dioxide. Note that the amount of CO increases at 3000 K initially and then decreases rapidly as CO is consumed by combustion for the fuel lean condition. In addition, since complete combustion of brown coal occurs for fuel lean conditions, the entropy gain due to the density change has a negligible effect on brown coal oxidation in this condition.  (Table 1); (e-h): S10 system (Table 1).
As shown in Figure 12, the variation of gaseous products with time is plotted. The number of oxygen molecules in S11/S12 decreases faster than that in S7/S8 as water molecules facilitate the reactions between oxygen molecules and brown coal. The CO/CO 2 production in S11 & S12 is similar to that in S7 & S8 because the oxygen molecules that are abundant in systems that complete combustions take place. However, as demonstrated in Figure 12c,g, the amount of carbon oxide increases as density is reduced by half. This increase in number of CO molecules can be attributed to the fact that CO is more stable at high temperature (3000 K), which agrees with previous results [22]. As discussed above, because brown coal combustion for fuel lean conditions is almost complete, the entropy change due to density difference has a negligible effect on final gaseous products formations.  (Table 1); (e-h): S12 system (Table 1).

Conclusions
In this manuscript, using a combination of ab initio molecular dynamics simulation (CPMD) and ReaxFF reactive molecular dynamics simulation, the interactions between water molecules and brown coal molecular model were investigated. Our CPMD simulation results demonstrate that hydrogen bonds formed between water and oxygen-containing polar functional groups contained in lignite, especially carboxyl, dominate the interaction between water and coal. Our simulation results not only confirm the results of previous quantum chemistry calculations [8,9], but also provide extra statistical information of the brown coal-water systems as system entropy is absent in those quantum chemistry calculations. The comparison between our CPMD simulations results and our ReaxFF-MD simulation results suggests that the ReaxFF reactive force field is not accurate enough to describe the diffusive behaviors of water molecules on the brown coal molecular model. Although the low-gradient method might correct the ReaxFF reactive force field at room temperature, it is still difficult to simulate temperature dependent diffusion of water on lignite, since the low-gradient method is based on the system configuration at a fixed low temperature [57].
Our ReaxFF-MD simulation results reveal the mechanisms in which the brown coal molecular model is oxidized for a wide temperature range. By varying the system density, it has been observed that the generation of gaseous products is sped up by density decrease due to the entropy gain. Hence, our simulation results as well as previous ReaxFF-MD simulations of coal combustion might not be able to describe the spontaneous combustion of coal piles as ReaxFF-MD simulations mainly focus on low density systems. Moreover, Zheng et al. [62] demonstrated that many reaction pathways can only be detected using ReaxFF-MD with a large molecular model of coal. Therefore, in order to study the low temperature oxidation of brown coal to elucidate the reaction mechanisms of coal spontaneous combustion, a more condensed phase of large brown coal molecular model with air is required.
The combination of CPMD and ReaxFF-MD simulations allows us to partially overcome the limitations of conventional force field, as we seek to study systematically the correlated water transportation and brown coal oxidation at various temperatures. How-ever, a modified version of ReaxFF reactive force field which can accurately describe the temperature dependent diffusion of gases (e.g., water) on lignite as well as a large brown coal molecular model with appropriate density are required to investigate the complicated reactions during coal spontaneous combustion/gasification. Author Contributions: Conceptualization, R.C. and G.W.; methodology, S.Y.; software, S.Y.; validation, R.C. and X.L.; formal analysis, S.Y.; writing, S.Y.; project administration, X.L. and G.W.; funding acquisition, X.M. and S.Y. All authors have read and agreed to the published version of the manuscript.